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Phonon scattering calculations predict the drag force acting on defects and dislocations rises linearly 
with temperature, in direct contradiction with molecular dynamics simulations that often finds 
the drag force to be independent of temperature. Using the Mori-Zwanzig projection technique, 
with no recourse to elasticity or scattering theories, we derive a general Langevin equation for a 
crystal defect, with full treatment of discreteness and non-linearity in the defect core. We obtain 
an analytical expression for the drag force that is evaluated in molecular statics and molecular 
dynamics, extracting the force on a defect directly from the inter-atomic forces. Our results show 
that a temperature independent drag force arises because vibrations in a discrete crystal are never 
independent of the defect motion, an implicit assumption in any phonon-based approach. This 
effect remains even when the Peierls barrier is effectively zero, invalidating qualitative explanations 
involving the radiation of phonons. We apply our methods to an interstitial defect in tungsten 
and solitons in the Frenkel-Kontorova model, finding very good agreement with trajectory-based 
estimations of the thermal drag force. 


Crystalline materials are invariably host to a huge 
population of defects such as dislocation lines, dislo¬ 
cation loops, vacancies, impurities and self-interstitial 
atoms. As is well known, the defect dynamics are 
typically a non-inertial mixture of drift and diffusion 
due to significant interaction of defects with thermal 
vibrations 1 . In this paper we describe in detail a recently 
reported method 2 for evaluating and understanding the 
interaction between thermal vibrations and crystal 
defects, resulting in a treatment of defect mobility that 
resolves an acute failing of phonon scattering theories. 

The dynamical law used to describe defect motion typ¬ 
ically balances deterministic forces from elastic interac¬ 
tions f against a viscous drag —7V, where v is the de¬ 
fect velocity (possibly that of a node on a dislocation 
line 3 ), whilst 7 is the drag coefficient, giving v = f/7. 
In dislocation dynamics literature 3 the drag coefficient 
7, sometimes labelled by B , is related to the mobility 
M = I/7, so that x = M f; in this paper we will use 7 
throughout. A common issue with this purely viscous dy¬ 
namical law is the absence of any temperature and thus 
thermal fluctuations. It has recently been shown 4,5 that 
to correctly capture the highly stochastic trajectories of 
nano-scale defects and kink-bearing dislocation lines seen 
in experiment 6,7 one must also add a stochastic thermal 
force ?y(f) to the defect equation of motion, giving the 
Langevin equation 

v = f/7 + 7?(t)/7, (1) 

where 77(f) is a white noise defined by the ensemble 
averages 8 

(77(f)) = 0, (r](t)r](t')) = 2k B T7(5(f - f'), (2) 

derived from the fluctuation-dissipation theorem 9 . The 
stochastic thermal force 77 is particularly relevant for 
nanoscale defects and small dislocation loops as they 


only respond to stress gradients 10 , meaning that the 
elastic force f is often negligible and the stochastic 
force dominates, giving the purely diffusive equa¬ 
tion of motion x = 77(f)/7 with a diffusion constant 
D = lim i _ >00 (a; 2 (t))/2f = k B T/7. For long dislocation 
lines the elastic force / typically dominates, giving an 
expected velocity of (x) = lim t ^. 00 (x(f))/f = f/7, as the 
expected stochastic force (77) = 0 vanishes. In either 
case, the value of 7 sets the time scale of defect motion 
and defines the dynamics of ensembles of interacting 
defects, and thus is a critical parameter in any disloca¬ 
tion mediated process such as post-irradiation annealing 
of defects, plastic deformation or the ductile-brittle 
transition. 

Whilst significant effort has gone into evaluating the 
deterministic elastic forces acting on crystal defects, 
there has been relatively little theoretical effort focused 
on the evaluation of the drag parameters 7 that play 
such a crucial role in any simulation of defect dynamics. 
All theoretical estimates use the results of Nabarro 11 , 
Eshelby 12 and others 13,14 , who calculated dissipation 
rates due to scattering by thermal phonons. In such 
treatments the drag force —yv is calculated to be 
proportional to the phonon density, meaning the drag 
coefficient is predicted to rise at least linearly with 
temperature in the classical limit, i.e. 7 ~ 7 w k B T. 
This so-called ‘phonon wind’ relationship is universally 
invoked, but has at best partial, qualitative, agreement 
with molecular dynamics (MD) simulation, agreeing with 
some simulations of dislocation lines in FCC metals 15,16 
but not identical simulations in BCC metals 17-21 , 
where the analysis of the average velocity under stress 
reveals 7 = f/(x) = 7o + 7wk B T, where 70 is indepen¬ 
dent of temperature. Most dramatically, simulations 
of nanoscale dislocation loops 4,22 , screw dislocation 
kinks 5,19 and highly mobile self-interstitial crowdion 
defects 23 in BCC metals, reveal a diffusivity that rises 



2 


linearly with temperature, meaning 7 = k B T /D = 70 
is independent of temperature. This surprising result 
is in complete disagreement with the phonon wind 
theory, which as we show below explicitly forbids the 
existence of a temperature independent component 
7o- Qualitative explanations for 70 claim dissipation 
arises due to phonon radiation as the defect passes over 
the Peierls migration barrier 14,20 lack any quantitative 
foundation and can clearly not account for the presence 
of 70 when the Peierls barrier is extremely small, as 
is the case for e.g. kinks on screw dislocations 24 and 
crowdions 25 in BCC metals. 

The lack of any quantitative theory for the presence of 
a temperature independent drag force highlights a need 
for a fundamental, quantitative understanding of the in¬ 
teraction of crystal defects with thermal lattice vibra¬ 
tions. In this paper we detail a recently reported 2 treat¬ 
ment of the defect drag force that fully accounts for the 
non-linear, discrete character of a real crystal and does 
not rely on elasticity and phonon scattering theories. We 
introduce defects as general localized deformations in the 
atomic configuration of a simulated crystal (section II), 
deriving a stochastic equation of defect motion using the 
established Mori-Zwanzig technique 26 (IIB). It is shown 
that in appropriate limits the drag coefficient 7 is equal 
to the integrated time autocorrelation of the defect force, 
divided by temperature, a Green-Kubo relation for the 
defect force. By taking averages over a quadratic Gibbs 
distribution of appropriately defined fluctuations defined 
in the subspace orthogonal to coordinates of the mov¬ 
ing defect itself, we derive an analytic expression for the 
autocorrelation of the defect force (III) and find that in 
general 

7 = 7 o+k B T 7 w , ( 3 ) 

where 70 and 7 W are constants. This central result pro¬ 
vides the a quantitative, theoretical demonstration of a 
temperature independent drag force for a general crys¬ 
tal defect with full treatment of non-linearity and dis¬ 
creteness, explaining the findings of many independent 
simulation studies 5,19,20,22,23 . 

We analyse the form of our explicit analytical expres¬ 
sion for 70 (section III), showing that the vibrations 
orthogonal to the defect motion are in general not 
vibrational modes of the crystal, by which we mean 
that vibrational displacements orthogonal to the defect 
motion produce a force that has a component parallel to 
the defect motion, even to linear order in the vibrational 
displacements. This means that phonons are not perfect 
oscillators perturbed by anharmonic couplings and 
consequently phonon momentum cannot be conserved 
in a defective crystal, invalidating the results of any 
phonon-defect scattering theory. These conclusions 
are explicitly tested on the discrete Frenkel-Kontorova 
model 2 ' (IV), a very special case as the continuum limit 
is an integrable system, with the ‘defect’ becoming a 
soliton 28 . In this special case 70 is predicted to vanish 


in the continuum limit, which we test in numerical 
simulations. Although the Peierls barrier is observed 
to vanish as we approach the continuum limit, we find 
that 70 remains, as fluctuations in a discrete lattice 
are always distinct from those in the continuum. This 
result shows that the presence of 70 is not related to any 
migration barrier effects, further invalidating a qualita¬ 
tive explanation of 70 based on phonon radiation 14,20 . 
In addition, the fact that 70 remains even in such an 
idealised system is a testament to the generality of our 
results. 

In the final section (V), we apply our approach to a 
realistic crystal defect, the 1/2(111) interstitial crowdion 
in Tungsten. Our approach offers a clear method to 
calculate the force acting on a defect directly from the 
inter-atomic forces, allowing a force autocorrelation and 
thus 7 to be calculated in molecular dynamics simulation 
runs and also by evaluating our analytical expression 
in molecular statics, finding very good agreement with 
typical measurements of 7 from diffusive simulation 
trajectories. It is hoped the present work will allow 
a better understanding of dissipation in crystalline 
systems and highlights the extent to which real crystal 
defects are fundamentally discrete and anharmonic 
objects that can never be assigned a conserved energy 
and momentum in a bath of canonical phonons. 


I. A QUALITATIVE SUMMARY OF 
CLASSICAL PHONON SCATTERING THEORY 

Previous attempts to calculate an effective drag force 
on crystal defects all aim to calculate the scattering 
cross-section that defects offer to thermal phonons. 
Whilst we will not go into a detailed derivation of these 
calculations 11,13,14,29 we provide a qualitative summary 
of how phonon scattering theory forbids a temperature 
independent drag parameter, i.e. 70 = 0, in conflict with 
extensive MD simulation results. 

A founding assumption of phonon scattering theory is 
that defects and phonons can each be considered ini¬ 
tially free particles, with quadratic and linear dispersion 
relations respectively. Drag occurs due to interactions 
through well defined scattering processes. In such a pic¬ 
ture, the drag parameter 7 = 70 + 7 W k B T + 72 (k B T) 2 +... 
represents a thermally averaged cross section to all scat¬ 
tering processes. As the phonon density is proportional 
to k B T in the classical limit, higher order terms in 7 will 
be the result of scattering processes involving a greater 
number of phonons. It turns out 11,14 that the terms 
of order (k B T) n involve n + 1 phonons, as illustrated 
in Figure 1 . For 70 the scattering process is therefore 
particularly simple- the absorption or emission of a 
single phonon. 

To see why a phonon scattering approach invariably 
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FIG. 1 . A cartoon of phonon scattering by a dislocation. 
Higher order scattering processes contribute terms of higher 
order in temperature to the effective drag parameter 7. As 
demonstrated in the text, phase space limitations result in the 
first order term 70 = 0 vanishing for subsonic defect speeds. 
The second diagram represents a typical scattering process 
which leads to the well-known 11,12,14 ‘Phonon Wind’ relation¬ 
ship 7 = 7wk B T. 


predicts that 70 = 0, consider a crystal with acous¬ 
tic wave speed c, lattice constant a and atomic mass 
m. Phonons have an energy E = \p\c and momentum 
\p\ < h/a <C toc, whilst dislocations have an energy 
E = P 2 / 2 fh and momentum P, where to > m /5 is 
the dislocations’ effective unit core cell mass. As we are 
far from any shock front and under typical mechanical 
stresses the dislocation speed is subsonic (far below the 
wave speed) i.e. P -C me. With a final dislocation mo¬ 
mentum of P' = P + p the energy balance thus reads 

P 2 (P + v ) 2 

—+ |p|c=——-,=> P = mc + p/ 2 ~mc, ( 4 ) 


which clearly violates the requirement that P <C me. As 
a result, the phase space available for such a one phonon 
scattering process vanishes and we are forced to conclude 
that 70 = 0, meaning 7 = yikpT + 72kpT + ... must rise 
at least linearly with temperature. This very reasonable 
argument has been invoked by all authors investigating 
the scattering of defects by thermal phonons 12,14,30 , but 
unfortunately is in conflict with many MD simulations 
of defect motion, as mentioned above. 

We have seen that the conclusion 70 = 0 is a direct 
consequence of assuming defects are canonical objects- 
initially freely moving particles that interact via scatter¬ 
ing processes, with a well defined and conserved energy 
and momentum. In the following sections we show that 
this founding assumption is always false as defect motion 
is fundamentally discrete and anharmonic. Also, the vi¬ 
brational modes of the crystal always change with defect 
position, meaning phonons cannot be defined as canoni¬ 
cal objects with conserved energy and momentum, inval¬ 
idating all the calculations of phonon scattering theory. 
Instead, to obtain a correct expression for the defect drag 
force we must treat the defect position and velocity sim¬ 
ply as functions of the full crystal configuration, deriving 
a general equation of motion using the established coarse 
graining techniques of the Mori-Zwanzig formalism 26 . 


II. DERIVING AN EQUATION OF DEFECT 
MOTION 

A. Defects in a discrete crystal 

The state of a classical crystal of N atoms is entirely 
defined by the N atomic positions x, £ R 3 and velocities 
x, £ R 3 , where i £ [ 1 , AT]. It will be convenient to repre¬ 
sent these positions and velocities by two 3 N-dimensional 
vectors X,X £ R 3N , constructed from the atomic posi¬ 
tions and velocities with the tensor sum 

X = (xi,x 2 , ...,x N ), X = (xi,x 2 , ...,x N ). (5) 

Molecular dynamics algorithms assign a potential energy 
V (X) for the system and then integrate Newton’s equa¬ 
tion ?uX = — W(X) under thermodynamic constraints 
to generate classical dynamics. It is clear that this 
inherently discrete system cannot exhibit the elastic 
singularities that form defects in a continuum system. 
Instead, defects are localized deformations which can 
be described by assigning a set of M < N nodes, with 
positions r £ R M and velocities r = v £ R M . We 
emphasize that r, v have no associated equation of 
motion, unlike the position and velocity of a particular 
atom. The dynamics of these nodes will be a projection 
of the Newtonian dynamics of the whole crystal. Many 
methods exist for determining the defect position and 
velocity, including analysis of the centrosymmetry 
parameter 3 , Voronoi analysis 31 and finding peaks in the 
time averaged potential energy 5 . In the following we 
do not constrain ourselves by the choice of a particular 
method for the identification of the defect and only re¬ 
quire that any method used is consistent and repeatable. 

The number of nodes one should assign to a crystal 
defect will obviously increase with the defect size. For 
example, to fully capture the configurational complexity 
of an extended dislocation line it is necessary to assign a 
node to each atomic plane normal to the dislocation line 
direction 5 . However, in the present work we focus, in 
the interest of clarity, on very simple defects which can 
be represented by a single node that only moves along 
a single direction, allowing us to set M = 1 , such that 
r = r £ R and v = v £ R. 

At zero temperature, there exists a well defined atomic 
configuration X=U(r) that minimizes the potential en¬ 
ergy of the crystal V (X) for each value of defect position 
r. These configurations are what is found in, for example, 
nudged elastic band calculations 32 . At a finite tempera¬ 
ture, these minimum energy configurations will be aug¬ 
mented by the displacements due to thermal vibrations 
$ £ R 3N , meaning the crystal configuration X at any 
given instant can be expressed as 

X = $ + U(r), X = $ + v< 9 r U, (6) 

where d r is the partial derivative with respect to r whilst 
keeping 3 > constant. For the 3 -/V-dimensional vector U(r), 
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U(r)-U(r)+6rdU(r) 





9U m+3 (r) 


FIG. 2. An illustration of the defect translation vector d r U 
for a localised ‘hump’ in a chain of ‘atoms’. The vector d r U 
describes the individual atomic displacements that correspond 
to an infinitessimal defect migration at zero temperature. 


which depends only on a single variable r, d r is equiva¬ 
lent to simple ordinary differentiation and can be readily 
evaluated. In the general case, we use the differential op¬ 
erator V = <9/<9X, to define <9 r through its action on a 
(possibly vectorial) function A(X) as 


d r A(X) = lim 

S —>o 


A(U(r + <$) + *)- A(U(r) + $) 


<9A 

dr 


= (9 r U • V) A, 


(7) 


which is clearly the infinitesimal change in A due to mo¬ 
tion along <9 r U. As illustrated in figure 2, the vector 
d r U is very important as it defines the directions of 3N 
atomic displacements necessary for defect migration and 
thus singles out the defect dynamics from the dynamics 
of the crystal as a whole. The mathematical treatment 
of d r U can be introduced using the trivial case of the de¬ 
fect being a physical atom. Assigning an index j to the 
atom, the 3N-dimensional vector 9 r U consistent with the 
variation of, say, the y-coordinate yj of atom j. is sim¬ 
ply d yj U = (0,0...1, 0...0) € R 3N . This vector picks out 
the relevant atomic coordinate from a full crystal config¬ 
uration. All the components of vector (0, 0...1,0...0) are 
zeros, apart from the one component in 3 j — 1 position, 
which is equal to 1. The set of the 3 N — 1 other direc¬ 
tions orthogonal to (0,0...1, 0...0) in this case are defined 
by the 3 N x 3 N matrix Diag(l, 1...0,1...1), where zero 
is in 3 j — 1 position. The action of this matrix on an 
arbitrary 3A-dimensional vector generates a vector that 
is orthogonal to (0,0...1, 0...0). 

For a localised deformation such as a crystal defect, the 
vector <9 r U now defines a direction in the 3A-dimensional 
space of atomic coordinates associated with the motion 
of a defect from r to r + Sr, as illustrated in figure 2. The 
matrix defining the space of all directions orthogonal to 
9 r U is now given by 


<9 r U <g> <9 r U 
<9 r U • <9 r U ' 


( 8 ) 


To illustrate the point, consider the action of this matrix 
on an arbitrarily chosen vector Y, 


Z = 



<9 r U ® <9 r U 


Y = Y 


<9 r U • Y 
<9 r U ■ <9 r U 


<9 r U. 


(9) 


Clearly, Z is orthogonal to <9 r U since 

a 'U Z = 9,U.Y-|Y|F(S rU .Y) = 0. (10) 

To see how this relates to the description of a crystal 
at finite temperature, we first note that the space of vari¬ 
ables including phonon displacements, velocities, and a 
single defect position and velocity $ ® © r © v has two 

more dimensions, i.e. 6 N + 2, than the 6A-dimensional 
space of atomic coordinates and velocities X © X. In or¬ 
der to ensure the same number of degrees of freedom in 
both coordinate sets we require the vibrational displace¬ 
ments $ to be orthogonal to the direction of displace¬ 
ments 9 r U caused by defect motion. This defines a 3N-1 
dimensional hypersurface of thermal displacements 33-35 

<9 r U • $ = 0. (11) 

The constraint (11) on the thermal vibrations $ can 
be derived by varying defect position r to minimize the 
quadratic deviation |X—U(r)| 2 , giving the minimum con¬ 
dition 9 r U ■ (X — U) = <9 r U • $ = 0. For a given value 
of defect position r, we can now vary v = r to mini¬ 
mize the quadratic deviation |X — v5 r U| 2 , resulting in 
<5|X — v<9 r U| 2 = (5v[<9 r U • (X — v<9 r U)] = 0, which in 
combination with the second of equations (6) gives 

<9 r U • $ = 0. (12) 

The constraints (11) and (12) emphasize that by as¬ 
signing a defect position and velocity to the crystal 
configuration, the 3N-dimensional vectors 3>, describ¬ 
ing the vibrational displacements of each atom in the 
system are restricted to lie on the 3N-1 dimensional 
hypersurfaces defined by (11) and (12) 36 . As each 
hypersurface will change as the defect moves (as <9 r U 
changes) the vibrational co-ordinates retain a depen¬ 
dence on r, although this dependence on r does not 
feature in the differential operator d r defined in (7) as $ 
is held constant. The representation of $ and 4> will be 
analysed in more detail when we calculate vibrational 
expectation values. 

To obtain a defect dynamical equation, we project 
the conventional atomic equations of motion mX = 
-VF(X) onto the vector 9 r U defining the direction of 
defect motion. Forming a scalar product of the equa¬ 
tions of motion with vector 9 r U, we obtain md r U • X = 
—(9 r U • W(X). We emphasize that as r is not a canon¬ 
ical co-ordinate we do not expect that the equation of 
motion for the defect to resemble a Hamiltonian equa¬ 
tion of motion; we aim to derive from the true equa¬ 
tions of motion mX = -VF(X) for the atoms an ex¬ 
pression that defines the dynamical law that defect po¬ 
sition r obeys. Differentiating the second of equations 
(6) with respect to time, and noting that r = v, we find 
that X = + v 2 <9 2 U + v9 r U. Exploiting the constraints 

(11),(12) we arrive at 

tov = -<9 r U • W(X) - (<9 r ro)v 2 /2 + v • • <f>, (13) 


<9 r U • <9 r U 
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where to = md r U • <9 r U is the effective mass of the defect 
and d T ifi = 2md^\J ■ <9 r U. Similar equations are known 
in other dynamical quasiparticle theories 28,34 . Following 
the established approximations, we neglect the ‘hydrody¬ 
namic’ term —v9 2 U • $ and the effective kinetic energy 
gradient — d r mv 2 /2. This is justified as we are only con¬ 
sidering defects with subsonic speeds |v| <C c and small 
migration barriers. For the effective kinetic energy gradi¬ 
ent, as the defect core is wide when the migration barrier 
is small (see below), the effective mass is of order 1/c 37 
and varies very little with position, meaning — 9 r TOv 2 /2 
is (at most) of order (v/c) 2 <C 1. The ‘hydrodynamic’ 
term —vd 2 U • <I> will have a vanishing expectation value 
as (<I>) = 0 (we will define this expectation value precisely 
below), which in turn means that this term will only ap¬ 
pear as second order or higher, which for a wide defect 
core will also give a contribution of order (v/c) 2 <C 1. 
Neglecting these terms gives a final dynamical equation 
for the defect coordinates of 


mv = -d r V = —<9 r U • W(X). (14) 


In appendix B we show that in the absence of the ther¬ 
mal vibrations and under a weak applied stress, the force 
on a crystal defect given by (14) is identical that fa¬ 
mously derived by Eshelby 38 . However, the main results 
of this paper concern the interaction of crystal defects 
with thermal vibrations; in order to perform the analyti¬ 
cal manipulation in section III we also require dynamical 
equations for the vibrational coordinates. In the same 
spirit as above, we analyse the full equation of motion in 
the subspace of all directions orthogonal to <9 r U, defined 
through the matrix (8), arriving at 


?n3> = — 


<9 r U 8) <9 r U 
(0 r U • <9 r U) 


• VF = -V*V. 


(15) 


where we have introduced the differential operator in the 
subspace of vibrations V$. As for d r in (7), V# can be 
defined through its action on a function A as 


V $ A 


0A\ 

d$) 


r 


d r U ® d r U 
(3 r U • <9 r U) 


• VA. 


(16) 


The two dynamical eciuations (14) and (15) clearly show 
a close parallel with the standard classical eciuations of 
motion toX = — W(X). However, the differential oper¬ 
ators (7) and (16) are defined through their relation to 
the direction of defect motion 9 r U in the 3N dimensional 
space of crystal configurations rather than the differenti¬ 
ation of a function with respect to a coordinate. Whilst 
we have now derived a general dynamical equation for 
the defect co-ordinates (14), they still retain an explicit 
dependence on the vibrational degrees of freedom. The 
next section uses the Mori-Zwanzig projection technique 
to replace the vibrational coordinates by a statistical dis¬ 
tribution, producing a closed but stochastic equation of 
motion for the defect. 


B. Removing the vibrational coordinates by the 
Mori-Zwanzig method 

From the form of the equation of motion for the defect 
(14) it is clear that the potential energy V(U(r) + <&) 
couples the evolution of the coordinate of the defect 
and vibrational coordinates. This is what is required 
for the frictional force to exist. The general dynamic 
relationship between these coordinates is therefore highly 
complex, but in the present work we only consider defects 
with low migration barriers moving at subsonic speeds. 
This allows us to assume that the defect coordinates 
may be considered as slowly varying compared to the 
vibrational coordinates, an approximation the validity 
of which we will explicitly prove later when we calculate 
the defect force autocorrelation in MD simulations. 

It is well known that defect migration barriers are 
directly related to the width of the defect core 39 . The 
wider is the defect core the lower is the migration 
barrier. This is because if a defect possesses a wide 
core, defect motion induces only small individual atomic 
displacements; in Peierls’ seminal paper 39 and many 
subsequent treatments 5,40 it has been shown that defect 
migration barriers decay exponentially fast with the 
defect core size. As here we consider highly mobile 
subsonic defects with very small migration barriers, we 
are therefore only concerned with wide defect cores, 
exhibiting broad, in comparison with the lattice pa¬ 
rameter, maxima in (<9 r U) 2 . We now use this fact to 
provide a heuristic argument justifying the assumed 
timescale separation. At a finite temperature, over 
a Debye period td ~ a/c ~ O.lps, where a is the 
lattice parameter, the displacements of any atom due 
to thermal vibrations have the oscillation amplitude of 
~ TD-y/ksT/m. Since the defect speed is approximately 
v ~ y/k^T/m <C c, the displacement of any one atom 
due to defect motion over a time interval td will be 
at most TDlldrUHooyffcBT/TO, where II^UHoo is the 
component of the greatest magnitude in 3 r U. These 
calculations imply that if II^UHoo <C |9 r U|, then the 
displacement due to defect motion will be much less 
than the magnitude of displacements due to thermal 
vibrations of atoms, which implies that the $ are 
effectively ergodic over a time-scale ~ tq where the 
defects essentially remains stationary. We note that the 
condition ||9 r U||oo ■C |<9 r U| amounts to a requirement 
that the deformation associated with the defect is spread 
over many atomic sites, which is always satisfied by 
highly mobile defects with a wide core. We therefore 
assume that vibrational displacements average to zero 
over periods of ~0.1ps whilst the displacements due to 
the defect structure are effectively static. Again, we will 
test this conclusion in MD simulation when calculating 
the defect force autocorrelation. 

This separation of time-scales can be exploited to inte¬ 
grate out thermal vibrations from the defect equation of 
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motion using the Mori-Zwanzig projection technique 26,41 . 
The idea is to derive a formal solution for the ‘fast’ co¬ 
ordinates <&, which may be substituted into the equa¬ 
tion for the ‘slow’ defect co-ordinates r, v. Crucially, by 
considering a distribution of initial conditions for the fast 
variables, we go from the micro canonical to canonical 
ensemble, introducing heat and stochastic fluctuations 
by only retaining statistical knowledge of the system. 
To actually do this we use a projection operator. The 
projection of some function A(r(t), v(t), <&(£), 3>(t)) is a 
conditional average over the fast variables weighted by 
some probability distribution function p; for systems with 
a well defined temperature (meaning the vibrations are 
to leading order harmonic 9 ) p is simply the conditional 
Gibbs distribution, given by 


p{ r, *) = exp (-0 [V (r, *) + |<J> • *]) /Z( r), 


(17) 

with the partial partition function Z( r) providing nor¬ 
malization. As defect velocity v only appears in the ki¬ 
netic energy mv 2 /2+TO$-4>/2, it does not couple directly 
to the vibrational coordinates and so does not appear in 
(17). As discussed above, although the vibrational dis¬ 
placements and $ appear as 3N-dimensional vectors, 
the constraints (11) and (12) restrict their allowed val¬ 
ues to a 3N-1 dimensional hypersurface orthogonal to 
(9 r U. Any conditional average over <E*, must therefore 
be taken on this hypersurface, resulting in a projection 
operator 


L = — U, which gives the time evolution of a general 
function A(r(t), v(£), 3>(f),4>(£)) of the crystal configura¬ 
tion through the relation gjA = LA. To derive the form 
of the Liouville operator we simply apply the chain rule- 

^A(r(t), v(£), $(£),#(£)) = v(t)d r A + v(t)d v A 

+ $•V$A + €>•V^A 

= LA, (21) 

which upon substituting in the equations of motion (14), 
(15) and requiring the identity to hold for any smooth 
function results in 

d r V 

L = — r ^d v + vd r -— • + $ • V*, (22) 

to in 

where we have introduced the differential operators V v 
and V^, the velocity space equivalents of <9 r and V#. 
Using Vx, which acts on the full crystal velocity X, these 
operators are defined as in (7) and (16), with d v = 3 r U • 
Vj and V* = (I - a r U ® <9 r U/(d r U • <9 r U)) • Vj. The 
formal solution of A A = LA reads 

at 

A(t) = exp ^(t — t')Lj A(t')i (23) 

with the special case A(t) = exp(fL)A(0). The exponen¬ 
tiated operator exp (tL) may thus be used to evolve any 
function, or projected function of the system coordinates 
in time. 


PA(t) = J A(r(i), v(t), 3>, &)p(r(t), 3>, 4>) 

x <5(<9 r U • $)<5(d r U • $)d$d$, (18) 

where the delta functions signify that we integrate over 
the hyper-surfaces defined by the constraints (11) and 
(12). The normalization condition on p is 

Z(r) = 

x (5(<9 r U • $)<5(<9 r U ■ $)d$d$. (19) 

This normalization condition eliminates any potential ar¬ 
bitrariness associated with the choice of the argument of 
delta functions in (18). To emphasize that the projection 
is a conditional expectation value, we will also employ the 
notation 

PA(t) = (A;r(t),v(t)). (20) 

This notation reflects the fact that P projects any 
function A(r(f), v(t), 3>(t),4>(t)) onto the space of func¬ 
tions of (r(f),v(£)); by definition, functions /(r(t),v(f)) 
only depending on (r(£),v(t)) are left unchanged, i.e. 
P/(r(£),v(£)) = /(r(£), v(£)). 

To obtain the key results as directly as possible it is 
expedient to use the anti-Hermitian Liouville operator 


The first step in our derivation is to partition the dy¬ 
namical evolution into a subspace of ‘fast’ and ‘slow’ vari¬ 
ables. We have already seen that the projection opera¬ 
tor P projects any function A(r(t), v(f), $(£), 4>(£)) onto 
the space of functions which only take (r(t),v(t)) as ar¬ 
guments. The projection operator is idempotent, i.e. 
P 2 = P, so a function /(r(t),v(t)) that only depends on 
(r(t), v(t)) is unchanged, i.e. / = Pf. We can also define 
the complimentary projection operator Q = I — P, which 
projects any function A(r(t), v(t), <!>(£), <&(£)) out of the 
space of functions which only take (r(t),v(f)) as argu¬ 
ments. It is simple to show that PQ = QP = P — P = 0. 
As a result, a function /(r(£),v(t)) only depending on 
(r(t),v(£)) vanishes under Q , i.e. Qf = QPf = 0. Just 
as P projected out defect coordinates, Q eliminates func¬ 
tions that do not have dependence on the vibrational 
coordinates. In this way, we can think of the ‘slow’ dy¬ 
namics being selected by P and the ‘fast’ dynamics being 
selected by Q. 

The two projection operators P and Q give a natural 
way to partition the defect equation of motion tov = 
-d T V(r(t), v(£), 3>(£), $(£)) = -d r V(t), using P + Q =1 
to write 

mv(t) = -(d T V-,r(t)) - Qd T V(t), (24) 

separating the force into a force expectation —Pd r V ( t) = 
— (d T V;r(t)) that depends only on the defect coordinate 
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and a remainder —Qd r V(t) = ( d r V;r(t )} — d T V(t) which 
clearly vanishes under P as PQ = 0. Using the Liouville 
operator L to extract the time evolution, the defect force 
then becomes 

mv(t) = —e t ^(d I V;r( 0)) — e t ^Qd r V( 0), (25) 

where exp(tL) is used to evolve the system from its con¬ 
figuration at a time t = 0. Following a standard method 
we now use an integral identity 26 

e tL = e^ 1 + f d se {t ~ s)L PLe s Q L , (26) 

Jo 

that is designed to arrive at the desired result as quickly 
as possible. This result can be checked by noting the in¬ 
tegrand is equal to e tL -^(e~ aL e s< ^ L ). Using this identity 
in (25) gives 

fnv(t) = —(d r V ; r(t)) + rj Q (t) + [ ds(Lr) Q (t);r(t - s)), 

Jo 

(27) 

where we have defined a ‘noise’ term 

VQ (t)=e^ L Qd r V( 0), (28) 

and a ‘memory’ term (Lr)Q[t)\?(t — s )) which by the def¬ 
inition of the projection operator as a conditional expec¬ 
tation may be written 

{L7] Q (t);T(t -s)} = J p(r(t - s), $, 4>)Lp Q it) 

x (5(«9 r U • $)<5(<9 r U • 4>)d$d$. (29) 

The noise term PQit) evolves under e t( ^ L , which can 
be identified as a ‘fast’ dynamics operator as functions 
/(r) only of the defect position r are not propa¬ 
gated: QLf(r ) = v<9 r / — Pvd r f = 0, meaning that 
exp(fQL)/(r) = /(r). 

Equation (27) is already reminiscent of a non- 
Markovian Langevin equation with a noise term that van¬ 
ishes under P (meaning it has an expected value of zero) 
and a memory term that is non-local in time. However, 
to complete the derivation for our purposes we need to 
put the integrand of the memory term into a more useful 
form. Using the anti-Hermitian property L = —U of L 
we can act on the distribution function p(t — s) in (29) 
instead of r]Q(t). We do this in stages for clarity. First, 
we recognise that 

d T Z(v) = -P(d I V;r)Z(r,v), (30) 

meaning that d r p(t ) = —pQd r V, which as V v p(t) = 0, 
gives the identity 

p(t — s)L = Dp(t — s) = — Lp(t — s) 

= pit — s)v(t — s) Qd r V{0). (31) 


Defining a second noise term 

p[t) = e tL Qd r V(0), (32) 

which evolves under dynamical operator e tL , unlike PQit) 
that evolves due to e t( ^ L , we can now re-express (29) as 

(LpQit)-, r(t - s)) = C(s; t)v(f - s), (33) 

C(s;f) = Pipit - s)p Q it);rit ~ s)), (34) 

which defines the memory kernel U(s;f) as the time cor¬ 
relation of the defect ‘noise’ forces p(t) and PQ{t). In 
general, the memory kernel retains a dependence on the 
absolute time t rather than the time difference s, as the 
kernel may depend on the defect position r(f — s) used in 
the expectation average (34). However, as we are consid¬ 
ering cases where the defect has small migration barriers 
we expect the average frictional force to vary little with 
the defect position, i.e. U(s;t) = C(s) an assumption 
that we have found to be validated in numerical simu¬ 
lations. We can now give a formally exact equation of 
motion for the defect coordinates as 

fnx.it) = ~id T V ; r(t)} + pcjit) — [ dsU(s)v(t - s). (35) 

Jo 

To evaluate the various terms in this equation it would be 
convenient to replace pQ(t) by pit), thereby eliminating 
the reduced evolution operator e t ® L present in ? 7 g(t) in 
favour of the full evolution operator e tL present in pit), 
as it is e tL that evolves the system in MD simulation and 
thus may be readily calculated. Such an approximation 
is possible when the defect dynamics occur on a slow time 
scale t over which the vibrational coordinates lose coher¬ 
ence, as in this limit the defect is effectively stationary 
over the dynamical range of interest, meaning e tL and 
e t Q L evolve the system in effectively the same manner. 
In this limit the rate of change of the defect variables 
iPL) is quantified by the small parameter r _1 , meaning 
in particular that exp (tQL) = exp(fL) + 0(r -1 ) and, due 
to the presence of Q in p and p q, 

CO) = Pipit - s)p(t)\ r it - s)) + 0(r -3 ) (36) 

is accurate to order r -3 , see Ref. 26 . This is a considerable 
simplification, because pit) = (d r V(f); r(f)} — <9 r U(t) is 
simply the fluctuating part of the defect force, meaning 
the memory kernel C(s) is simply the autocorrelation of 
the fluctuating force. As stated above, in this timescale- 
separated regime we expect C(s) to decay rapidly with s. 
As a result we take the Markovian approximation, where 
the noise force is replaced by a delta correlated white 
noise and v(t — s) is taken out of the integral in (35). The 
time integration of C(s) is formally extended to infinity 26 
though in practice we shall see that a lower limit, usually 
of order r, is more appropriate for computational and 
physical reasons 42 . With these approximations we finally 
obtain a Langevin equation for a crystal defect that reads 

mv(t) = -(<9 r U;r(f)) -7 v(t) + p{t), (37) 



though it is typical (and usually legitimate 5 ) in disloca¬ 
tion dynamics to neglect the inertial term rhv(t), which is 
valid when the potential energy landscape is slowly vary¬ 
ing over the thermal length ^keT/my 5 ’ 8 . In this limit 
we finally obtain the overdamped Langevin equation 

v(i) = -(< 9 r V;r(f )}/7 + j?(f)/ 7 , (38) 

which is identical, after identifying —(d T V;r(t)) with the 
elastic forces, to (1) given at the start of this paper. In 
equations (37) and (38), {r){i)fj{t')) = 2 keT 7 < 5 (t — t') is a 
purely Markovian white noise force 8 and 7 is the time in¬ 
tegral of the defect force autocorrelation divided by tem¬ 
perature, viz. 

/*oo 

7 = /? / C(s)ds 

Jo 

/*oo 

= /? / P(p(t-s)p(t);i(t-s)). (39) 

Jo 

Equation (39) is our main result, a Green-Kubo type 
relation 9 that equates the defect force autocorrelation to 
the defect drag parameter 7 . We have thus achieved a 
central aim of this paper, to produce a Langevin equation 
for a crystal defect starting from the classical equation 
of motion for the 3N atoms of a defective crystal. To see 
the connection of (39) to typical Green-Kubo relations, 
consider the Fokker-Planck equation for a free particle 
with position x and momentum p: 

d t p(x,p, t ) = ( p/m)d x p - 7 pd p p + k B Tjdpp. (40) 

The last term has the form of a diffusion operator in 
momentum space, with a ‘momentum space diffusivity’ 
keTy 8 . Therefore, just as the well-known expression 
D = f Q (v(s)v(0))ds relates the velocity autocorrelation 
(v(s)v(O)) to the real space diffusivity D, equation 
(39) relates a force autocorrelation / 0 °° C(s)ds to the 
momentum space diffusivity keTy. 


The real utility of (39) is that we can calculate the 
memory kernel directly from simulated trajectories of 
the stochastic defect force —d T V(t) = — <9 r U • W. In 
the present case of small migration barriers the expected 
force — (d T V;r(t)) will depend only on long range, slowly 
varying elastic forces which will be uncorrelated with 
thermal vibrations, contributing a term solely dependent 
on the defect position. This term will not contribute to 
pit) = — Qd r V(t ) as it will be removed by Q , meaning 
that we can invoke a standard ergodicity assumption to 
write 
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force autocorrelation, we now evaluate it in two ways. 
In section III an analytic derivation of the vibrational 
expectation values gives the expected autocorrelation 
(39) as function of curvatures of the potential energy, 
which can be evaluated in static calculations in section 
IV for the Frenkel-Kontorova model and in section V 
for realistic crystal defects. This approach gives detailed 
insight into the form of the interaction between crystal 
defects and thermal vibrations, explaining the existence 
of a temperature independent drag coefficient, but 
numerical evaluation typically requires a full eigenvector 
decomposition which is expensive for large systems. 
In section V we present a much more efficient method 
that directly evaluates the defect force — <9 r U • W in 
molecular dynamics simulations, extracting (41) through 
analysis of the defect force time series. 


III. ANALYTIC DERIVATION OF 7 

To derive an expression for 7 we must remove 
by calculating the vibrational expectation values (45), 
integrating over the hypersurfaces <9 r <I> = 0 and <9, <I> = 0 
given in (11) and (12). In particular, we want to calculate 
the the defect force autocorrelation {p(t — s)p{t)) for use 
in (39), where 77 (f) = —d r V(t) + (d r V) is the fluctuating 
component of the defect force. To do this, we first expand 
both the potential energy and the defect force in powers 
of around X = U(r). The expanded potential energy 
reads 

V = Vo(r) + -<&-VJ,Vo ■ 3> + T^y3>■ V 4 Vo • 3?• <& + ..., (42) 

where the subscript Vo indicates that all the partial 
derivatives are calculated at X = U(r). The first vibra¬ 
tional derivative V 3 .V 0 = 0 as by definition (<&) = 0 at 
zero temperature. 

In order to perform an integration over we must 
parametrise the hypersurface of allowed values, namely 
all directions orthogonal to 3 r U. To do this we define 
a basis set of 3 N — 1 orthonormal vectors {efc(r)}™j _1 
which are all orthogonal to 9 r U, i.e. 9 r U • e*,(r) = 0 V 
k. In order for this condition to hold for all values of 
the defect position the {e^} will vary with r, though as 
we perform the integrals for a given value of the defect 
position this will not affect the result of our integration. 
With a set of 3 N — 1 vibrational amplitudes 
and momenta ' we are free to express $ and 

as 

3.Y -i 3JV-1 

*« = V ^ ^ Pk (r) 7 (43) 

k =1 k =1 


Having obtained the Green-Kubo type relation (39) which explicitly shows the dependence of the vibrational 
for the defect drag parameter in terms of the defect co-ordinates on r. We can use this parametrization to 
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evaluate (18) through a change of variables in the inte¬ 
gration. As the {e^.} are by definition orthonormal the 
Jacobian in (18) is simply unity for <1? and l/m 3 ^' 1 for 
< 1 ?, giving 

r 3N-1 

PA(t)= / A(r,v, {<f>k,Pk})p d<j)k(dp k /m), (44) 

k 

= (A;r(t),v(t)). (45) 

It is convenient to also define the differential operator 


This choice of {e^,} as the eigenbasis of V\V means that 
Vl^Vo = 8 qr vuu 2 is diagonal. If there is no external 
potential the crystal will have global translational sym¬ 
metry, meaning three eigenvectors will always be closely 
related to a global rigid translation with extremely small 
eigenvalues; we will omit these when calculating vibra¬ 
tional expectation values. After this choosing the vi¬ 
brational basis vectors to be the eigenvectors of the <I>- 
Hessian V\V the expanded potential energy takes the 
simpler form 


= e k ■ V, (46) 

which, as the {e*,} span all directions orthogonal to <9 r U 
gives the identity 


eh 2 

V = V 0 (r) + -f muJ q + 


3! 


r *Vo 


(53) 


whilst the defect force expansion reads 


3N-1 

I — d r U 0 d r U/(d r U ■ d r U) = e k ® e fc . (47) 

fc =i 

As the differential operator V# is by construction or¬ 
thogonal to <9 r U (see equation (16)) we can therefore ex¬ 
pand as 

3N—1 

v* = ^2 <t, k - (48) 

fc=l 

In the remaining manipulations we will employ the Ein¬ 
stein summation convention 43 , where repeated indices 
are taken to be summed from 1 to 3N — 1. For exam¬ 
ple, the expansion of differential operator can now 
be written 


V#=e^V 0 fc . (49) 

With these identities we can now rewrite the expanded 
potential energy, using the Einstein summation conven¬ 
tion, as 


V = l / o(r) + 




V 0 + .., (50) 


where again all the partial derivatives are evaluated at 
X = U(r). This last expansion immediately suggests an 
expedient choice for the basis vectors {efc}: the eigen¬ 
basis of the matrix of second derivatives in (42), the ‘<f>- 
Hessian’ 


LJ 0 d r U\ 

u-a r u )' 

(51) 

With knowledge of V(X) and <9 r U this real, symmetric 
matrix is readily constructed and can be fully diagonal¬ 
ized to obtain 3N orthonormal eigenvectors and eigenval¬ 
ues. The <f>-Hessian V|,V is distinct from the full Hessian 
V 2 H as the vector of defect motion 3 r U is by construc¬ 
tion an eigenvector with an eigenvalue of identically zero, 
meaning that the 3 N — 1 remaining eigenvectors form an 
orthonormal set {efc(r)}, orthogonal to 9 r U, with ‘eigen- 
frequencies’ {wfc} such that 

• e fc = (Vl k V 0 )e k = mule k . (52) 


VlV = 


I- 


<9 r U 0 d r U 

au • au 


VT' I- 


ai 


ay = ay 0 +a 9 av^y 0 + ^av^v 0r y o +.., (54) 

where in both cases the Einstein summation convention is 
used. Practically, the mixed derivatives in (54) are eval¬ 
uated by contracting the full tensorial derivative V"U, 
evaluated at X = U(r), once with aU and n — 1 times 
with the {efc}. It is important to note that although 
V 09 y = 0, there is no requirement that the mixed deriva¬ 
tives vanish; moreover, these terms give the coupling be¬ 
tween the defect and thermal vibrations. To proceed, we 
now truncate V to quadratic vibrations in the Gibbs dis¬ 
tribution (17). Our choice {efc} means that the Gibbs 
distribution is a product of Gaussians, namely 

p(r, {^>fc,Pfc}) = e -^[Vo(r)+rr,E,«^;/2+E,p a l /H/Z(r) > 

(55) 

meaning expectation values may be calculated using 
standard identities such as Wick’s theorem 43 . This trun¬ 
cation clearly neglects any thermal expansion arising 
from the purely vibrational anharmonicities V(}- 3 y. In 
a previous publication 2 we systematically included these 
terms to produce an expression for 7 up to linear order 
in temperature. It was shown that our main result, the 
anomalous temperature independent drag coefficient 70 
in equation (59) below, is unaffected by these additional 
terms. 

To analytically evaluate the defect force autocorrela¬ 
tion and hence the memory kernel C(s), we evolve the 
vibrational coordinates from a fixed r. This is justified by 
the time-scale separation, as C(s) is expected to decay to 
incoherence before r changes appreciably, an assumption 
we will test in MD simulation. To do this we truncate the 
vibrational force to linear order in (15), differentiate (43) 
and invoke the timescale separation to neglect terms of 
order v, obtaining an approximate vibrational equation 
of motion of 


<iq = (56) 

We emphasize that this equation of motion only holds 
over short time periods where the defect co-ordinate r 



10 


is effectively stationary and the defect force autocor¬ 
relation C(s) is expected to decay to incoherence, an 
assumption we will test explicitly in section V. Over 
longer periods of time r, and hence the vibrational basis 
{efc(r)}, will change; nevertheless, on the short timescale 
of present interest (56) may be solved with initial con¬ 
ditions ((/>g(0)(/> r (0)) = (k B T /mu 2 )S qr , (^(0)^(0)) = 0, 
allowing us to explicitly evaluate time dependent expec¬ 
tation values of {(f>k} in terms of the two point correlation 
functions 43 


{<t> q {t)(t> r(0)) = - -S qr cos(oj q t). (57) 

mojq 

As appropriate for non-conservative dynamics, the corre¬ 
lation function is evaluated using only initial conditions 
and consequently is closely related to the retarded Green’s 
function 44 

G q {t) = ®(t)P{<l> q {t)<l> q (Q)) = cos(aV), (58) 

where 0(t) is the Heaviside step function. Finally, us¬ 
ing Wick’s theorem 43 to simplify high order expectation 
values {(j) q (j) r <f) s ), (<j> q <j>r<f>s(f>t) in the force autocorrelation, 
we perform the elementary Gaussian integrations to give 
our main result 
/*00 

7 = / (d r V^V) 2 Gq(t)dt 

Jo 

7 rji nOQ 

+ ^J Q d r Wl^VG r (t)G s m^Vdt 

+ ^ f~ Wt'VGqWGriWVl^Vdt (59) 

We see that the friction coefficient takes the form 7 = 
7o + k B Ty w , with a temperature independent component 


7o 



(drVt'V) 2 

mcjg 


cos(oj q t)dt. 


(60) 


The temperature independent drag coefficient 70 couples 
the defect to thermal vibrations via a mixed quadratic 
derivative 


drV<t, k V = drU- V 2 V-e k , (61) 

where V 2 V is the full Hessian. This is the first time, 
to our knowledge, that a theoretical derivation has pro¬ 
duced an expression for the widely observed temperature 
independent drag force acting on a crystal defect. 

The existence of this coupling means that vibrational 
displacements $ that are orthogonal to defect motion 
(as 3 r U • $ = 0) can still induce a force on the defect 
(through —3 r U • V 2 K • $) to linear order in the vi¬ 
brational displacements. As a result any theory that 
assumes that thermal vibrations are harmonic oscillators 
perturbed by anharmonic terms will never capture 70 


as it is assumed a priori not to exist. In our approach, 
although we approximate the dynamics of 3> by (14) 
over short timescales, our general derivation made no 
assumption on the form of the crystal potential energy 
and thus allows any form of coupling to exist. 

The temperature dependent drag coefficient k B T 7 w 
couples through the mixed cubic and quartic deriva¬ 
tives <9 r V 2l F and d T V^V. As typical scattering theo¬ 
ries only include cubic anharmonic 13,14,45 only the sec¬ 
ond term in (59) is directly comparable. In a rough 
approximation, if we consider the Hessian as approxi¬ 
mately V 2 V ~ 22 tow 2 cubic anharmonicities become 
V 3 G ~ 2 ^ 2 m.w(Vw), meaning the cubic coupling term 
is approximately (V 3 K) 2 /mw 2 ~ (22(^ w )/ w ) 2 ; roughly 
comparable to square of the Griineisen parameter 46 70 = 
20 fp-/ w, which measures the change of vibrational fre¬ 
quencies with hydrostatic pressure. As a result we con¬ 
clude that the temperature dependent frictional coupling 
is of order k B Ty w ~ k B T( 7 c ) 2 in agreement with other 
approximate analyses 45 . The precise nature of the anhar¬ 
monic coupling will be the subject of a future publication. 
However, in the remainder of our analysis we concentrate 
on the nature of the anomalous temperature independent 
7o- 


IV. A SPECIAL TEST CASE: KINKS IN THE 
FRENKEL-KONTOROVA MODEL 

We have seen in section I how phonon scattering 
theories forbid a temperature independent drag param¬ 
eter, 70 = 0 , as thermal vibrations are taken to be 
the phonon modes of the perfect lattice with, before 
the inclusion of anharmonic vibrations, a conserved 
energy and momentum. To better understand how the 
assumptions of phonon scattering theory contrast with 
the present approach, we investigate a rare system which 
supports localised deformations but also admits, in the 
continuum limit, a full analytic evaluation of all the 
terms (<9 r U, {e g }, V^V, etc.) used in our analysis. This 
will highlight the precise effect of canonical thermal 
vibrations on the motion of crystal defects. 

The Frenkel-Kontorova model 47,48 is built from a set 
of N harmonically coupled nodes with one-dimensional 
positions X = (X 0 , X 2 ,..., Xjq-i) € R w , sitting in a sinu¬ 
soidal ‘lattice’ potential of period b. The potential energy 
reads 27,48 


VFK(X) = ^.| (Xitl ~f l ~‘‘ ) 

(62) 

where a is the equilibrium spacing of the harmonic cou¬ 
pling and k, v are energy densities that control the rela¬ 
tive strength of the coupling and ‘lattice’ potential. The 
kinetic energy is simply 22 , a(p/2 )(A,) 2 , where ap is the 
node mass and the system is typically completed with 


av sin 


X,; 



11 


periodic boundary conditions Xjv = Xo + Na. Dynamics 
are classical, i.e. produced by 


a/DQ = —V Xi Vfk(X). 


(63) 


The continuum limit takes the discrete chain of nodes 
{X„} to an elastic line X(x), such that X„ —> na+X(na). 
As a — y 0, (X.j_^i — X ? ; — of / a — y VX(x), Net — y L — y oo 
and o, — S ► J da. As with all such limits, directions in a 
vector space now become functions in a Hilbert space 43 , 
with inner products becoming integrals. In this limit we 
obtain the so called sine-Gordon model 


V[X] = £ dx^(VX(x)j 2 


+ v sin 7r 


X(x) 


(64) 


with a kinetic energy f Q dx(/x/2)(X(x)) 2 . It is well 
known that with boundary conditions X(x + L) = X(a) + 
a the sine-Gordon model support solitons of the form 
X(a) = U(a,r), where the defect configuration U(a,r) 
reads 5 


U(a, r) = - + — arctan ( sinh ( 

2 7T V V 


x — r 

w 


(65) 


The soliton width w is given by w = ( b/2Tr)^/2n/v) 
(we assume L w) and the soliton energy is E = 
(2bITt)\J2 kv, which is clearly independent of solition po¬ 
sition, i.e. the migration barrier is identically zero. It 
is a simple matter to give the direction of defect motion 
(now in Hilbert space) as 


c> r U(a, r) = ——secli ( --- 

TTW V W 


giving an effective mass of 

rL 


m = 


y f (d r U(a, r)) 2 da = 
Jo 


2d 2 y 


( 66 ) 


(67) 


The sine-Gordon model is an integrable system, which is 
particularly remarkable in that one can derive exact ana¬ 
lytical forms for the set e 9 (a,r) of the vibrational modes 
orthogonal to 9 r U which remains valid for all values of 
r. The vibrational modes are given by 48 


e g( x > r ) = 


/ 7 iqw 
\ L 


i tanh 


x — r 


))exp( 


7 xqx 

i— - lUJ n t 

Li 


( 68 ) 

where ui q = \Jrtly\J 1 + ( nqw/L ) 2 , meaning that ther¬ 
mal vibrations can be included as 


$(x,r)= / dq(f q e q (x,r), 

Jo 


m4>(a,r) = / dqp q e q (x,i), 
Jo 


(69) 

(70) 


thereby giving analytic expressions for all the quantities 
used in the theory developed in section III. It is simple 


to show that the orthogonality conditions (11), (12) are 
satisfied as 48 

f da<9 r U(a,r)e 9 (a,r) = 0. (71) 

Jo 

However, 9 r U and the {e 9 (a,r)} are also eigenmodes of 
the full Hessian V^^a') = 5(x — a')V 2 H(a), where 

tt 2 v sinh 2 ((a — r)/rc) 1 

\7 2 V(x) = kV 2 + -^- ) - L - L -. 72) 

cosh 2 ((a; — r )/w) 

3 r U has an eigenvalue of identically zero 

[ S7 2 V(x,x')d, U(a',r)da' = V 2 H(a)<9 r U(a,r), 

Jo 

A) = °’ ( 73 ) 

7 TW° / 


2 sinh 2 ( Yff ) — cosh 2 ( Yff ) / a 


2 / x— r 


cosh 3 (^) 


TTW TTW 


whilst the vibrational modes e 9 (a,r) have eigenvalues of 
/iw 2 , namely 

[ X7 2 V(x, x')e q (x', r)da' = V 2 H (x)e q (x, r) 

Jo 

= yw 2 e q {x, r). (74) 

We thus see that the sine-Gordon model is a very special 
case of the general derivation of section III, where we 
found that thermal vibrations $ could be decomposed 
in the eigenmodes |e g (r)} of the ^-Hessian (51) evalu¬ 
ated at a given value of r. In the special case of the 
sine-Gordon model, whilst it is still true that the ther¬ 
mal vibrations $(a, r) can be decomposed into the eigen¬ 
modes e g (a,r) of the $-Hessian V %V{x,x') 1 in this case 
the 4>-Hessian V|H(a,a') is equal to the full Hessian 
X7 2 V(x,x'). To see this, we first construct the Hilbert 
space ‘matrix’ Q(x,x', r) of all directions orthogonal to 
(9 r U, namely 

rM / ^ su ^ d Y \J(x, r)3 r U(a; / , r) 

Q{x,x ,r) = S(x — x) — L - -j—, (75) 

Jo (<9 r U(a,r)) da- 
meaning that the <f>-Hessian V %V(x,x') reads 

nL pL 

VlV(x,x')= / Q(x,y,v)V 2 V(y,z)Q(z,x’,r)dydz, 
Jo Jo 

(76) 

which is the continuum equivalent of the discrete <f>- 
Hessian (51). However, because <9 r U is an eigennrode of 
the full Hessian with eigenvalue zero, as shown in (73), 
we have that 

[ Q(x : y,i-)V 2 V(y,z)dy = \7 2 V(x,z), (77) 

Jo 

and similarly for the other inner product, giving the final 
result that 

V|V(a, x') = 5{x - x')V 2 V(x) = V 2 V(x, a'), (78) 
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showing the equivalence of \/%V(x,x') and V 2 V(x,x') 
for this system. In this rare case of an integrable system 
we thus have a convergence of the approach developed 
here and the assumptions of phonon scattering theory, 
indeed we can use the correct expression for the vibra¬ 
tional modes (rather than plane waves of scattering the¬ 
ory) and treat the vibrational momentum and energy as 
conserved quantities (rather than, as in our approach, 
simply an expression of the crystal configuration con¬ 
tingent with a given r). As a result, in this case the 
eigenmodes e q (x,r) of the full Hessian are orthogonal to 
the defect motion mode i9 r U, as shown in equation (71). 
From our analytical result (60) we have seen that the 
temperature-independent drag parameter 70, is a sum 
over all the modes of the <f>-Hessian V|K weighted by a 
coupling term (61) <9 r U • V 2 V • e/ c . In the Sine-Gordon 
model, it is simple to show that this coupling vanishes, 
viz. 

<9 r U • V 2 y • e k -> [ dxd r U(z,r)v2He 9 (a;,r), (79) 
Jo 

= u ;/w(.,.Ww)=o, (so) 

which in turn means that the temperature independent 
drag parameter 70 = 0 vanishes. This result is a direct 
consequence of the integrability of the sine-Gordon 
model, a rare property that allows one to identify a 
general eigenmode expansion encompassing the direction 
of defect motion <9 r U and all other orthogonal vibrational 
modes e q (x,r), valid for all values of the defect position. 
As the vibrational modes are now modes of the full 
Hessian rather than just the <I>-Hessian, the amplitude 
(j> q of each vibrational modes can be considered as an 
isolated harmonic oscillator, independent of any defect 
before the inclusion of anharmonic vibrational terms, 
with a conserved energy and momentum. 

Phonon scattering calculations also assume that 
phonons are harmonic oscillators with a conserved 
energy and momentum, interacting with defects only 
after the inclusion of anharmonic vibrational terms. But 
we have just seen that for these assumptions to hold 
in our general approach we are forced to consider only 
integrable systems, rare non-linear models that typically 
only exist in one dimension 4 '. It is therefore very 
interesting that only in this somewhat abstract limit 
does our general approach come into agreement with the 
central tenet of any phonon scattering calculation. This 
highlights the ‘artificial integrability’ that scattering 
calculations must impose on a problem in order to have 
a well defined phase space for phonon scattering. In 
general we do not expect these assumptions to hold and 
thus we anticipate 70 / 0 to be a general feature of 
defects in any realistic crystalline system. 

These considerations promote the Frenkel-Kontorova 
(62) model as an interesting test case to calculate 
70, as by better approximating the continuum limit 


we reach the Sine-Gordon model, where by equation 
(80) 70 should vanish. To this end we have performed 
static calculations on the Frenkel-Kontorova model 
(62), which in appropriate regions of parameter space 
also supports soliton-like ‘kink’ deformations under 
boundary conditions X„ + jv = X„ + b. To investigate the 
continuum limit we scale the line tension n and lattice 
barrier v in inverse proportion, i.e. (k, v) — \ (cm, v/a) 
such that the kink energy ~ is constant and the 

kink width ~ \J k/v — > a.yj k/v increases; in this way the 
deformations of the chain become smoother and better 
approximate the continuum field of the sine-Gordon 
model, whilst as the kink energy is roughly constant we 
can compare results taken for different values of the kink 
width. 


It is well known 5,4 ' that for kink widths greater than a 
few node spacings (i.e. w/a> 2 the system state is very 
well approximated by the discretised soliton kink solution 
X„ = U„(r) = na + U (na — r), where U(na — r) is given 
by (65). However, when we substitute this solution into 
the form of the total energy (62) the discrete summation 
retains a dependence on kink position, giving a system 
energy of approximately 5 


E( r) = (2b / it)\/2kv + -—-—exp(— nw/a) sin 2 (7rr/a), 

a 

(81) 


being the continuum result (2b/n)y/2ni' plus a periodic 
term, of period a, that decays exponentially fast with the 
kink width w. The presence of this periodic migration 
barrier is an accepted signature of discreteness effects 
in crystalline systems, with the exponential decay with 
defect width famously derived by Peierls 39 to explain 
the low critical shear stress of dislocations. 


The kink migration barrier decays rapidly with 
increasing kink width as a broader kink profile causes 
the summation in the total energy (62) to better 
approximate the integral (64) of the continuum limit. 
As this integral has no dependence on kink position 
the kink migration barrier vanishes. It is interesting 
to ask whether we expect a similar dependence for the 
temperature independent drag parameter 70. 

Whilst the kink migration barrier arises from the kink 
formation energy being a function of defect position in a 
static system, the expression (60) for 70, in contrast, is 
a sum across all the fluctuations in the system, weighted 
by a coupling term (<9 r U • V 2 V ■ e q ) 2 . 

In the continuum limit of the Sine-Gordon model 
we have seen (equation (80)) that this coupling 
term vanishes, as the {e 9 (a;,r)} are eigenmodes 
of the full Hessian V 2 V(x)6(x — x'), namely 
that fy V 2 V(x,x')e q (x',r) = fj,co 2 e q (x,r). In 
this case the discrete summation that approxi¬ 
mates f^V 2 V(x, x')e q (x',r) is the inner product 
V 2 K • e q , where the summation is over all elements 
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[(e g )i, (e g ) 2 , •■•(e q )jv] of the vector e q . However, unlike 
the kink profile, the fluctuations e q in a discrete system 
have no requirement to vary slowly between lattice sites, 
regardless of the kink width or any defects the same sys¬ 
tem. As a result we cannot say that an increasing kink 
width suppresses discreteness effects in the fluctuations, 
meaning we expect y 0 to decay much more slowly than 
the migration barrier as we attempt to approximate the 
continuum limit. 

We emphasize that this is a general feature; even if 
discreteness effects are negligible for the migration of 
stable defect structures that possess a broad displace¬ 
ment field, the fluctuations around a crystal defect of 
any size will always be sensitive to discreteness as they 
have no requirement vary slowly across lattice sites and 
thus approximate a continuum solution. 

To test these conclusions we have evaluated 70 for 
a kink in the Frenkel-Kontorova model, employing the 
techniques detailed in section V below. We construct 
and fully diagonalize V%V to obtain {e q ,uj q } and U(r), 
allowing us to evaluate the integrand in (60), which 
is then integrated with respect to time (see below) to 
produce 70 . In agreement with our earlier assumption 
that the frictional force is essentially independent of 
defect position, the actual value of 70 was found to vary 
very little with the precise position of the defect r. 

The results of these calculations are shown in Figure 
3. We see that whilst the temperature independent drag 
coefficient does indeed eventually vanish as the contin¬ 
uum limit is reached (with increasing kink width) it does 
so at a much slower rate than the kink migration bar¬ 
rier, which rapidly becomes indistinguishable from zero. 
We therefore conclude that a temperature independent 
friction parameter is robust, arising due to the coupling 
of a localised defect to fluctuations in a discrete system. 
These fluctuations give an important discreteness effect 
even in finely interpolated approximations to integrable 
systems, where they are expected to vanish, and typ¬ 
ical static signatures of discreteness such as migration 
barriers are indistinguishable from zero. This demon¬ 
strates that a temperature independent drag coefficient 
70 arises even when migration barriers vanish, invalidat¬ 
ing a proposed ‘radiative damping’ mechanism that ar¬ 
gues 70 arises due to a defect radiating phonons when 
traversing migration barriers. As 70 is present even in 
this idealised system it is reasonable to expect this phe¬ 
nomenon to arise for all crystal defects, which we inves¬ 
tigate in the next section, applying the developed theory 
to realistic crystalline systems. 


V. NUMERICAL EVALUATION OF 7 

The previous analytically tractable test case allowed a 
detailed analysis of our main result (60), but the present 



FIG. 3. Migration barrier E m i g and temperature independent 
friction coefficient 70 as a function of kink width w, relative 
to their values at w = a. We see the migration barrier decays 
exponentially, whilst the temperature independent friction co¬ 
efficient remains almost unchanged even for kink widths much 
larger than the lattice parameter a. This highlights the fact 
that 70 is a quantity owing its existence to the discreteness 
of the system, which is present even when the energy cost of 
discreteness is low. (Colour Online) 


work was motivated by the observation that the drag co¬ 
efficient 7 ~ 70 is independent of temperature for many 
nanoscale defects in atomistic molecular dynamics (MD) 
simulations. For such nanoscale defects the drag coeffi¬ 
cient is typically extracted under zero applied stress by 
calculating a diffusion coefficient 4,5,23 from the defect tra¬ 
jectory. It is well known that for a particle of mass m 
at temperatures above any migration barriers present we 
have the Ornstein - Uhlenbeck relation 8 


2 


k B T, 


1 — exp ( — — 


(-U 

\ m / 


(82) 


that shows the transition from ballistic behaviour 
(x 2 (t))/t 2 ~ mk B T at short times to diffusive behaviour 
at long times with the Einstein relation 9 


D = 


lim 

t—t OO 


2 1 


k B T 

7 


(83) 


Equation (83) will be the benchmark value for 7 against 
which we test our results. In this paper, we have simu¬ 
lated a 1/2(111) interstitial crowdion defect in Tungsten, 
which is known to posses an anomalous temperature in¬ 
dependent drag coefficient 7 = 7 o 22 - The simulations 
were performed with the LAMMPS package 49 using an 
EAM potential 50 for W by Marinica et al. 51 . For the cal¬ 
culation of D from (83) we used supercells of ~ 50,000 
atoms (30 x 30 x 30 unit cells) which were relaxed, ther- 
malized using a Langevin thermostat, and then run for 
l- 2 ns over a range of temperatures to extract a diffusion 
constant D, as shown in figure 4. As D was observed 
to rise linearly with temperature, we performed a least 
squares fit and used (83) to find a value of 

7 = 6.0 (7) eV • fs/A 2 , 


(84) 
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FIG. 4. The diffusivity of an interstitial crowdion defect in 
W at various temperatures. A linear fit gives D = kBT/ 70 . 
Inset: typical data for the mean square displacement against 
time. We see ballistic behaviour at very short times, leading 
to diffusive behaviour after ~ 0.5 ps.(Colour Online) 

where units are a reflection of the fact that 7 repre¬ 
sents the impulse density of a heat bath. We also note 
that with an effective mass rh ~ m/7, we can predict a 
‘ballistic-diffusive transition’ time-scale of 

Tfl 

— ~ 0.43 (5) ps, (85) 

7 

which is entirely consistent with the ( x 2 {t )) data in figure 
4, showing that the overdamped Langevin equation (38) 
can adequately capture the stochastic dynamics of crystal 
defects on timescales larger than a picosecond. 

A. Evaluation of C(s) in Molecular Dynamics 

We have developed a numerical scheme to determine 
d r U in dynamical MD simulation, allowing us to calcu¬ 
late the defect force 

-d T V{t) = —<9 r U • VV(i) (86) 

which can then be used to calculate the force autocorre¬ 
lation C(s). From the analysis of section IIB we derived 
the result 

/»00 

7 = / C(s)ds, (87) 

Jo 

thereby allowing us to calculate 7 from the trajectory of 
the defect force. The scheme is as follows. We first use 
the stable crystal defect structure at zero temperature 
to define the defect configuration U(r), using the nudged 
elastic band (NEB) barrier climbing technique 32 to 
generate U(r) for all values of the defect position within 


a lattice period parallel to the defect motion. All other 
lattice periodic images may then be generated simply 
by rigid translation. For the simple interstitial defect 
considered here we have found that to an extremely 
good approximation the set of NEB configurations may 
be generated from a single ground state configuration, 
which can offer computational efficiencies in the dynam¬ 
ical implementation of the method. 

As detailed in appendix A, as the deviation of U(r) 
from the configuration of a perfect lattice varies slowly 
across atomic sites, we can equate d r U(r) to the atom¬ 
istic strain field along the direction of defect motion. 
We can then use linear interpolation to evaluate U(r) 
and d r U(r) for all values of the defect position r within 
a lattice period parallel to the defect motion. The 
application of this method to extended defects will be 
presented in detail elsewhere, though we have also found 
good results for certain glissile edge dislocation lines 2 
and we expect this interpolation technique to be valid 
for edge-type defects with wide cores, as the atomic 
displacements are in the direction of defect motion and 
vary slowly with defect position. 

Once U(r) and d r U have been calculated, we run 
an ensemble of finite temperature, zero-stress MD 
simulations of the same system extracting the 3N- 
dimensional crystal configuration X(i), velocity X and 
force — W(X(t)) under micro-canonical conditions 5 . 
Whilst in principle we can extract this data every MD 
timestep (here lfs), we have found that data can be 
taken every 5 — 10 MD timesteps, still capturing the 
fastest vibrational modes in the force autocorrelation 
(~ 50fs) whilst gaining some efficiency from time 
coarse-graining the subsequent data analysis. For each 
extracted configuration we find, to within a small 
tolerance, the zero temperature configuration U(r) that 
minimizes the quadratic weight |(X — U(r)) • d r U| 2 . 

The choice of the quadratic weight |(X — U(r)) • d r U| 2 
is important, as its minimum should ideally vanish, as 
|(X — U(r)) • d r U| 2 = |$ • d r U| 2 = 0, regardless of the 
defect structure or temperature. We find that this can 
be achieved to a very good approximation at finite tem¬ 
perature, as shown in the inset of figure 5; the minimum 
value of the quadratic weight |[X — d r U] • d r U| 2 is typi¬ 
cally ~ 10 7 , meaning a threshold value of |(X — U(r)) • 
d r U| 2 < 10 -5 gives satisfactory accuracy in the defect 
position. Figure 5 demonstrates that after minimising 
| [X — U(r)] • d r U| 2 as described above, any crystal config¬ 
uration X(f) can be split into a defect configuration U(r) 
and a featureless, fluctuating field of thermal vibrations 
3?. This is an important justification of the claim (6) 
made at the beginning of section II that such a decom¬ 
position was possible at finite temperatures. 

We have found that the determination of U (r) at each 
time step can be significantly accelerated by using the 









15 



Atomic position along [111] (Units of b=V3a/2) 



FIG. 6. The defect position, velocity and force acting on the 
defect, extracted for an interstitial crowdion in W at T=300K. 
To our knowledge this is the first time a defect velocity and 
force have been extracted directly from the velocity and force 
vectors of a simulated crystal. (Colour Online) 


FIG. 5. Determination of U(r) from MD simulation of a 
1/2(111) crowdion in W at T=300K. Below: The deviation in 
1/2[111] bond length for U(r) and X(t). Inset: illustration of 
a 1/2(111) crowdion from 52 . Above: The thermal vibration 
vector $ = X(t) — U(r), which fluctuates around zero with no 
peaks, as expected. Inset: Logarithmic plot of the quadratic 
weight |(X — U(r)) • d r U| 2 for various values of r. We see a 
quadratic minimum of ~ 10“' which may be readily detected. 
(Colour Online) 

identity md r U • X(i) = mv(f) (see figure 6 ) to approx¬ 
imate r(t + St) by r (t) + Stv(t). We have found this 
effective preconditioning of the minimisation condition 
means a satisfactory minimum may be found extremely 
quickly, with our method becoming comparable in 
speed to the typical method of extracting a defect 
position through analysis of the atomic disregistry. The 
efficacy of this method is a confirmation of our earlier 
assumption that the defect position typically varies very 
little over a thermal vibration time scale of ~ O.lps 
~ 100 MD timesteps. 

As shown in figure 6 , the defect position trajectory 
obtained in this manner is typically much smoother than 
those found by analysing peaks in the potential energy 
or centrosymmetry parameter, as we fit a rigid defect 
profile to the whole configuration, removing fluctuations 
from the fitting procedure. However, the real benefit of 
our approach is that we may use 9 r U to project out the 
defect velocity v = <9 r U • X(f)/(<9 r U • 9 r U) and the force 
— d r V = —9 r U • W(X(t)) acting on a defect at each 
timestep. To the best of our knowledge this is the first 
time such quantities have been directly extracted from 


the atomic forces and velocities in an MD simulation. 

The defect velocity is important as it can be used to 
aid the fitting procedure as described above, but it also 
gives important data to test the assumptions in our ap¬ 
proach. We have calculated the trajectory average (v 2 ) 
over a variety of temperatures and found that it gives the 
equipartition value trike T to within 10%. Furthermore, 
the root mean velocity is typically 1 or 2 Burgers vectors 
b per ps, where b = \/2>a/2 is the translational period of 
the defect, meaning that over a vibrational period ~ 0.1 
ps the defect center moves at most 10 % of a migration pe¬ 
riod, supporting the timescale separation argument that 
we made at the start of section IIB. The defect force is, 
as expected, a rapidly fluctuating function of time with 
zero mean. 

With a defect force trajectory it is now a simple mat¬ 
ter to calculate the defect force autocorrelation, which is 
then divided by kgT to produce the memory kernel C{s). 
From equation (39), the time integral of C(s) should be 
equal to the defect drag parameter 7 . As shown in figure 
7, the initial peak of C(s) is essentially identical across a 
wide range of temperatures, with the subsequent signal 
oscillating around zero. Importantly, we see that C(s) 
decays to zero after ~ 0.05ps~ td/ 2 over which time 
the defect is observed (figure 6 ) to be essentially station¬ 
ary, again validating the assumption of timescale separa¬ 
tion we made at the start of section III. As the signal 
varies significantly and incoherently across different sim¬ 
ulations, and typically flattens as we increase the simula¬ 
tion time and system size, we limit the time integration 
to the first zero in the memory kernel, finding a value 
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FIG. 7. The defect force autocorrelation calculated from 
molecular dynamics simulations at various temperatures. The 
small supercell used for one of the T=150K autocorrelations 
contained ~ 3, 000 atoms, whilst the other data was taken 
from supercells containing ~ 10, 000 atoms. We see the ini¬ 
tial peak of the force autocorrelation divided by knT is es¬ 
sentially independent of temperature, giving an estimate for 
a temperature independent drag parameter 7 = 70 — 5.9eV- 
fs/A that compares well with the value of 70 = 6 .0(7) eV- 
fs/A obtained from measurement of the diffusion constant 
D = knT/ 7 o (inset). (Colour Online) 


for 7 = 70 that is in good agreement with 7 extracted 
from analysis of D = keT/y from long diffusive trajec¬ 
tories generated using larger simulation supercells, typi¬ 
cally containing 50, 000 atoms, as compared with 10,000 
atoms used for producing the force autocorrelation data. 
It is interesting to note that the peak of the force au¬ 
tocorrelation does not significantly vary across systems 
sizes and even for low temperatures (T=50K) where the 
very small lattice migration barrier ~ O.OleV ~ lOOfcg 
causes the defect not to migrate only very rarely over a 
simulation time of a nanosecond. The precise nature of 
these finite size effects and their role in data analysis is 
a subject of a separate study. 


Hessian 

U<g>d r U\ 

U • 9 r U ) ' 

(89) 

As this expression only involves the calculation of second 
order derivatives it can be readily evaluated from a zero 
temperature configuration. After minimising the system 
to obtain U(r), it is a simple matter to construct the 
matrix 


ViU = I - - 


<9 r U ( 8 > d r U 


au • au 


• V 2 U- I- 


<9,1 


dr 


drU (8) <9 r U 
<9 r U ■ dr U 


(90) 


which spans the 3 N — 1 directions orthogonal to <9 r U(r). 
We can then evaluate the Hessian matrix V 2 U of sec¬ 
ond derivatives (see appendix C) at X = U(r) to con¬ 
struct the ^-Hessian V|.U, defined in equation (51). 
Using standard LAPACK routines, we can then obtain 
the 3 N — 1 eigenvectors and eigenvalues { e q , 
which may be used to evaluate ( 88 ). The result of these 
calculations are shown in Fig. 8 . We see that the initial 
peak of the calculated force autocorrelation is in good 
agreement with the dynamical measurements and varies 
little when using a tiny system of 433 atoms ( 6 x 6 x 6 
supercell) compared to a larger system of 3457 atoms 
(12 x 12 x 12 supercell). The slight additional amplitude 
of the force autocorrelation for the larger cell can be ex¬ 
plained by analysing the coupling (<9 r U ■ V 2 V • e g )/mw 2 
of the defect to each vibrational mode, shown in the in¬ 
set of figure 8 . We see a remarkably similar profile in 
the number of modes that couple strongly to the defect, 
but that a larger system clearly has a greater number 
of modes that couple only weakly, which will all con¬ 
tribute a small amount to the summation in (60). The 
study of this vibrational coupling for a variety of de¬ 
fects will be the topic of a future work, but we finish 
this section by noting that our central result, the expres¬ 
sion ( 88 ) for the temperature independent memory ker¬ 
nel Cq(s), giving 70 = / Co(s)ds, closely resembles the 
famous Kac-Zwanzig memory kernel Cj<z(s) for a par¬ 
ticle connected to a bath of harmonic oscillators. The 
Kac-Zwanzig memory kernel reads 26 


B. Evaluation of C(s) in Molecular Statics 

In section III we derived an analytical expression (60) 
for the temperature independent 70 in the form 70 = 
/ 0 °° Co(s)ds, where the temperature independent mem¬ 
ory kernel Co(s) read 

CoO) = (d rU V J e i) cos ( w s ). (gg) 

muf. 

q <? 

The summation in ( 88 ) is over the set {e 9 , of 

vibrational eigenvectors and eigenfrequencies of the <&- 


Ckz(s) = cos (u g s), (91) 

q 9 

which is identical to Co(s) except that here the coupling 
constant A is a constant as opposed to (<9 r U • V 2 U • e q ), 
and the normal modes are typically assumed to be that of 
a continuous, homogeneous three-dimensional medium, 
i.e. ui = ck , with ~ w 2 of oscillators of frequency w, 
meaning the memory kernel (91) becomes precisely a 
delta function. It is interesting that crystal defects have, 
in realistic systems, a very heterogeneous coupling to 
their environment, but nevertheless the total coupling 
strength varies only a little as we consider larger systems. 
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FIG. 8 . The defect force autocorrelation calculated from 
molecular dynamics simulations at various temperatures. We 
see that the autocorrelation is independent of temperature, 
giving an estimate for a temperature independent drag pa¬ 
rameter 7 = 70 that can accurately predict the diffusion con¬ 
stant D = kaT/ 7 o (inset). (Colour Online) 

VI. CONCLUSIONS 

In this paper we have developed a general method to 
define, construct and evaluate the stochastic equation 
of motion for a crystal defect. By defining a crystal 
defect only as a localised deformation of a non-linear, 
discrete crystal of N atoms, we were able to identify a 
3 N dimensional vector 9 r U that gave the direction of 
defect motion in the configurational space of the crystal, 
i.e. the set of N atomic displacements that effect defect 
migration from a given defect position. This vector is 
readily evaluated in silico and may be used to extract a 
defect position, velocity, and force directly from the 3 N 
dimensional atomic positions, velocities and forces used 
in any MD simulation. This allows one to analyse much 
more than simply the defect trajectory with time and 
gives real fundamental insight into the defect dynamics. 

The vector 9 r U allowed us to define the set of all 3iV- 
1 mutually orthogonal directions as a subspace of ther¬ 
mal vibrations consistent with a given defect position, 
with an equation of motion for all of these quantities 
emerging from the atomic equations of motion. Using 
the Mori-Zwanzig technique, we then derived, under well 
defined approximations, a Green-Kubo relation (39) for 
the defect drag parameter 7 that equates 7 to the time 
integral of the defect force autocorrelation divided by 
temperature. A numerical implementation found that 
this method agreed well with traditional estimates of the 


diffusivity using the position trajectory, but that high 
quality data on the force autocorrelation could be ob¬ 
tained, due to the much smaller correlation time, from 
much smaller simulation supercells, offering an interest¬ 
ing method to extract long time behaviour from a rela¬ 
tively short dynamical simulation. 

Under the assumption of a well defined system temper¬ 
ature (Gaussian Gibbs distribution) we also derived an 
analytical result (59) for 7 in terms of the system poten¬ 
tial energy. We found that in general 7 = 70 + kBT 7 w , 
where 70 is a temperature independent term that is for¬ 
bidden to exist in all previous theoretical treatments 
based on phonon scattering formalism, but has been 
widely observed in MD simulations. We found that 70 
arises because the vibrations orthogonal to defect motion 
are in general not eigenmodes of the crystal. This was 
shown to violate a founding assumption of phonon scat¬ 
tering calculations, namely that thermal vibrations may 
be considered as the normal modes of a perfect crystal, 
which then interact with crystal defects through the in¬ 
troduction of anharmonic terms in the potential energy. 

In a real crystal, the presence of a defect drastically 
changes the normal mode structure, and in general this 
new set of normal modes will not be orthogonal to the 
direction of defect motion, meaning that the vibrations 
e q will couple to the defect even to quadratic order, i.e. 
through the quadratic coupling term e q ■ V 2 U • 9 r U 7 ^ 0. 

To test this conclusion we investigated an extreme case, 
the integrable sine-Gordon system, where analytical re¬ 
sults predict that in the presence of a soliton the nor¬ 
mal modes are all orthogonal to the direction of defect 
motion, meaning in turn that for this integrable model 
70 = 0. However, when simulating the sine-Gordon sys¬ 
tem, even with a very fine discretisation, we saw that 
whilst traditional static discreteness effects such as the 
Peierls barrier are effectively zero, 70 / 0 was still 
present, as fluctuations persist up to the smallest length 
scales in any discrete system, meaning they cannot be 
approximated away. As a result we expect the presence 
of 70 to be a general feature in the stochastic dynamics 
of any crystal defect. 

Future work will concern detailing the application of 
this approach to extended defects such as dislocation 
lines and glissile stacking faults, and investigating how 
these techniques can be applied to develop a stochastic 
equation of motion for defects that posses a large migra¬ 
tion barrier. 
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Appendix A: Calculation of <9 r U from a single atomic 
configuration 

In the present work we focus on the case where de¬ 
fects posses very small migration barrier, which invari¬ 
ably means that the atomic displacements caused by the 
defect vary slowly along the direction of defect motion, 
as a small migration barrier implies no individual atom 
migrates a large distance under defect migration 5,39 . We 
now exploit this slow variation to show how one can cal¬ 
culate <9 r U(r) from a single configuration U(r). Let the 
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defect move in a direction parallel to some lattice vector 
direction a, typically the Burgers vector b. Consider the 
positions of an individual atom u,(r) G M 3 , i G [1, N] and 
a neighbouring atom which we assign an index n(i ) G 
[1, TV] such that far from the defect u„(j)(r) — u, (r) —> a. 
We recover the full 3N-dimensional configuration U(r) 
through tensor addition, i.e. U(r) = ©jUj(r). By virtue 
of the translational symmetry of the crystal lattice we 
know that the displacement of atom n(i) for a defect po¬ 
sition r is equivalent to the displacement of atom i for 
a defect position r — a , namely u n o\(r) = u,(r — a) + a. 
For a slowly varying displacement field we may now ex¬ 
pand Uj(r — a) to give U;(r — a) ~ u,(r) — a9 r u, (r); com¬ 
bining these manipulations with the definition <9 r U(r) = 
©j9 r Uj(r) we finally obtain 

<9 r U(r) ~ 0 (u ra( j)(r) - u 4 (r) - a) /a. (Al) 


This result, valid when the atomic displacements vary 
slowly along the direction of migration, allows us to 
calculate <9 r U(r) simply through calculating a finite 
difference derivative for each atom and its nearest 
neighbour along a. 


Appendix B: Equivalence of — d r Vo to Eshelby’s 
configurational force 

As in the treatment of the main text, we consider a 
crystal at zero temperature containing a defect which 
may be described by a single position parameter r, such 
that the atomic configuration X = U(r) G R 3N . In the 
presence of weak external tractions T G R 3N which are 
assumed to be non-zero only far from the defect core, one 
can show that the defect force is precisely the configura¬ 
tional force on a crystal defect first derived by Eshelby 38 . 
As the applied tractions are weak, the total potential en¬ 


ergy may be written as 

y 0 (r) + T-U(r), (Bl) 

where the subscript Vo indicates that X = U(r). We can 
then be vary (Bl) with respect to r to obtain 

—<9 r U 0 = T • 9 r U(r). (B2) 

The requirement that T is weak and applied far from 
the defect core is to ensure the linearity in T of the to¬ 
tal energy (Bl); if the perturbation was non-linear the 
continuously parametrized set of minimum energy con¬ 
figurations U(r) would not be a global minimum in the 
presence of an external traction. We note that equiv¬ 
alence between elasticity and a fully non-linear discrete 
treatment is only expected to apply in this regime. To 
explicitly apply a surface traction, let T represent a force 
of ±Ah ■ a G R 3 per atom for two bounding planes E-t 
far from the defect core, where A is the area per atom. 
The defect force is now 

-d r V 0 = A E n • cr • <9 r Ui(r), (B3) 

Ui6S + ,E_ 

where 9 r Uj(r) G R 3 , ! G [1,1V] are the individual compo¬ 
nents of c> r U(r) for each of the N atoms in the system. 
To take a continuum limit we let the set of displacements 
u, (r) — x 3 , where x 3 is the nearest perfect lattice posi¬ 
tion, be interpolated by the continuum vector field of 
displacements u(x; r) to give 

— d r Vo = I dS'n • cr • <9 r u(a:), (B4) 

AeE + ,s_ 

in direct agreement with Eshelby’s result. By virtue of 
global rotational symmetry we can always ensure a is 
symmetric under permutation of indices, meaning we can 
always replace <9 r u(;r) in (B4) by its symmetric compo¬ 
nent, the strain field e = d r u(x) + (<9 r u(a;)) T . Under 
mild assumptions on the nature of the displacement field 
one may also continue Eslrclby’s derivation to obtain the 
Peach-Koehler dislocation force. 


Appendix C: Tensorial derivatives of an embedded atom potential 


The potential energy for a set of atoms {x 1 } interacting through an embedded atom potential are typically of the 
form 


U({x ! }) = U,[ A] - UyM 


= £ i E^') =E e 


c 




(Cl) 


where r lJ = |x* — x J | > 0 is the Euclidean distance between atoms i and j, A is a pair potential term and co, the 
keystone of the embedded atom method, represents the electronic density. In practice these potential terms are 
neglected once r lJ exceeds some cut-off radius r max . 
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1. Derivatives of the pairing function <j> 


The embedded atom potential (Cl) is built from pair-potential functions <^(|xj — Xj|) = <p 13 between pairs ij. To 
simplify later notation, we now define the first and second derivatives, which by the translational invariance of the 
argument will have permutation symmetry in the cartesian directions a, (3, 7 , e £ ( x , y , z)- 


Xa = 


W d(j) ij 


dxi 


= y M 

Aq 7 


do? 


W = M 
a/3 


dx 3 a 

dx% d X % 


dxi 


dx{ 


= -h 


dxi, 


(ij) 

(a/3)’ 




9y 


wi 


<9xjg 


= (5i* - 


(ij) 

(<*P) 


(C2) 

(C3) 

(C4) 


where [], () indicate the antisymmetric and symmetric permutation symmetry. Practically, the Cartesian derivatives 
X, ’h of <fi(r) are evaluated in spherical polar coordinates though we omit these standard results as there is quite 
enough algebra already. To aid the following, we also define the ‘reduced’ quantities 


Q C- n = 


Cl 


(C-n )! 


0 


C—n 


= ^2xc 


$ 


a/3 


= E** 


etc. 


(C5) 


dx^ 




'Y 'ik 
1 a(3 7 


etc. 


2. First and Second Derivatives 


With these definitions, we can immediately write 


dU c 

dxi 




de m 

dxi. 


C E 0 m _1 E W + S mn ) = CY, X^ (Bf- 1 + 9m" 1 ), 


(C 6 ) 


-fFh = C E %E( 0 E X + 0 - _1 ) + w -!) E xr ] (^ 0 E 2 + ©w 

WW m m dx J 0 dx P 

=ce $ ?^ - ^)( 0 f _i + 0 ™ _i ) 


+ C(C -1) £ xE ro] xjr (<% + <W©f “ 2 + xL^x?" 1 W + -we 


,\C —2 


Which in our reduced notation reads 


dUc = y* ©J- 1 y im 0 G_1 
dx i Xa KJ + / . ta '-h , 


+ 0 E X )+ExrxrW 2 


E xl'xffi- 2 + A, (wef- 1 + E W 0 ^ 1 + xi4 0 f“ 2 ,) 

')£F Ua,i8) \ m / 


(pa,q(3)£F(ia,j /3) 


(C7) 


(C 8 ) 


where the index-coordinate permutation sum is explicitly 

(pd,q/3) = (ia,j/3),(j/3,ia). (C9) 

As any analytic partial derivative must be invariant to the order of differentiation, it will be advantageous to group 
terms in higher order derivatives as sums over such index-coordinate permutations. The overbars on the Greek 
coordinate symbols are used to distinguish them from the coordinate symbols outside of such sums. 



